/*
 * kernels.cu
 *
 *  Created on: 2011-2-18
 *      Author: zhang
 */
#include <arguments.h>
#include <crs.h>
#include <stdlib.h>
#include <stdio.h>
#include <curand_kernel.h>
#include <pthread.h>
#include <sys/time.h>
#include <ell.h>

#define BLOCK_SIZE 128
#define BLOCKS 64
#define THREADS (BLOCKS*BLOCK_SIZE)

template<class T>
struct SharedMemory
{
    __device__ inline operator       T*()
    {
        extern __shared__ int __smem[];
        return (T*)__smem;
    }

    __device__ inline operator const T*() const
    {
        extern __shared__ int __smem[];
        return (T*)__smem;
    }
};
// specialize for double to avoid unaligned memory
// access compile errors
template<>
struct SharedMemory<double>
{
    __device__ inline operator  double*()
    {
        extern __shared__ double __smem_d[];
        return (double*)__smem_d;
    }

    __device__ inline operator const double*() const
    {
        extern __shared__ double __smem_d[];
        return (double*)__smem_d;
    }
};


__global__ void InitializeRand(curandStateXORWOW* state, long seed){
  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
  curand_init(seed,tid_global,0,&state[tid_global]);
}

/***********************
 * performs vector reduction of sum within a block
 **********************/
template<class T, unsigned int blockSize>
__inline__ __device__ void sumInBlock(T* sdata, int tid){
	  if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
	  if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
	  if (blockSize >= 128) { if (tid <  64) { sdata[tid] += sdata[tid +  64]; } __syncthreads(); }

	  if (tid < 32)
	  {
	      // now that we are using warp-synchronous programming (below)
	      // we need to declare our shared memory volatile so that the compiler
	      // doesn't reorder stores to it and induce incorrect behavior.
	      volatile T* smem = sdata;
	      if (blockSize >=  64) { smem[tid] += smem[tid + 32];}
	      if (blockSize >=  32) { smem[tid] += smem[tid + 16];}
	      if (blockSize >=  16) { smem[tid] += smem[tid +  8];}
	      if (blockSize >=   8) { smem[tid] += smem[tid +  4];}
	      if (blockSize >=   4) { smem[tid] += smem[tid +  2];}
	      if (blockSize >=   2) { smem[tid] += smem[tid +  1];}
	  }
}


template<typename T>
__inline__ __device__ T cukpmGetNormalRand(curandStateXORWOW_t* state);

template<>
__inline__ __device__ float cukpmGetNormalRand<float>(curandStateXORWOW_t *state){
	return curand_normal(state);
}

template<>
__inline__ __device__ double cukpmGetNormalRand<double>(curandStateXORWOW_t *state){
	return curand_normal_double(state);
}

/*********************
 * firstly perform generate vector R1 = H * R,
 * then perform v-v multiplication R*R1
 */

template<class T,unsigned int blockSize>
__global__ void cukpmCreateRandR1(curandStateXORWOW_t *state,T *r,T *r1,int size,T *mu)
{
	  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
	  int tid = threadIdx.x;
	  T mySum=0;
	  T *sdata = SharedMemory<T>();

	  for(int i=tid_global;i<size;i+=THREADS){
			  r1[i] = r[i] = cukpmGetNormalRand<T>(&state[tid_global]);
			  mySum += r1[i] * r[i];
	  }
	  sdata[threadIdx.x] = mySum;

	  __syncthreads();

	  //make summation within a block
	  sumInBlock<T, blockSize>(sdata, threadIdx.x);

	  // write result for this block to global mem
	  if (tid == 0)
		  mu[blockIdx.x] += sdata[0];
}

template<class T, unsigned int blockSize>
__global__ void cukpmCreateR2(ELL<T> *H,T* r1,T* r2,T* r,T* mu){
	int tid_global = blockDim.x * blockIdx.x + threadIdx.x;

	T mySum=0.0;
	T *sdata = SharedMemory<T>();
	for(unsigned int  k=tid_global;k<H->rowsNum;k+=THREADS)		// k-th row of matrix H
	{
		T product = 0.0;
		for(unsigned int j=0;j<H->colsNum;j++)
		{
			//H->elements[1]= 0;
			product +=  H->elements[j * H->rowsNum + k] * r1[(H->cols[j * H->rowsNum + k])];
		}
		r2[k] = product;
		mySum += product * r[k];
	}
	sdata[threadIdx.x] = mySum;

	__syncthreads();

	//make summation within a block
	sumInBlock<T, blockSize>(sdata, threadIdx.x);

  // write result for this block to global mem
	if (threadIdx.x == 0)
		mu[blockIdx.x] += sdata[0];
}

template<class T, unsigned int blockSize>
__global__ void cukpmCreateRn(ELL<T>* H,T* r1,T* r2,T* r,T* mu){

	int tid_global = blockDim.x * blockIdx.x + threadIdx.x;
	int tid = threadIdx.x;

	T mySum=0.0;
	T *sdata = SharedMemory<T>();
	for(unsigned int  k=tid_global; k<H->rowsNum; k+=THREADS)		// k-th row of matrix H
	{
		T product = 0.0;
		for(unsigned int j=0;j<H->colsNum;j++)
		{
			product += H->elements[j * H->rowsNum + k] * r2[(H->cols[j * H->rowsNum + k])] * 2;
		}

/*		product += H->elements[0 * H->rowsNum + k] * r2[(H->cols[0 * H->rowsNum + k])] * 2;
		product += H->elements[1 * H->rowsNum + k] * r2[(H->cols[1 * H->rowsNum + k])] * 2;
		product += H->elements[2 * H->rowsNum + k] * r2[(H->cols[2 * H->rowsNum + k])] * 2;
		product += H->elements[3 * H->rowsNum + k] * r2[(H->cols[3 * H->rowsNum + k])] * 2;
		product += H->elements[4 * H->rowsNum + k] * r2[(H->cols[4 * H->rowsNum + k])] * 2;
		product += H->elements[5 * H->rowsNum + k] * r2[(H->cols[5 * H->rowsNum + k])] * 2;
		product += H->elements[6 * H->rowsNum + k] * r2[(H->cols[6 * H->rowsNum + k])] * 2;
*/

		r1[k] = product - r1[k];
		mySum += r1[k] * r[k];
	}
	sdata[threadIdx.x] = mySum;

	__syncthreads();

	//make summation within a block
	sumInBlock<T, blockSize>(sdata, threadIdx.x);

  // write result for this block to global mem
	if (tid == 0)
		mu[blockIdx.x] += sdata[0];
}

template<class T, unsigned int blockSize>
__global__ void cukpmSumMu(T* d_mu,T* d_U, unsigned int m, T d){
	int tid = threadIdx.x;
	T *sdata = SharedMemory<T>();

	for(int b=blockIdx.x;b<m;b+=gridDim.x){
		sdata[threadIdx.x] = d_mu[b * blockDim.x + threadIdx.x];

	__syncthreads();

	//make summation within a block
	sumInBlock<T, blockSize>(sdata, threadIdx.x);

	  // write result for this block to global mem
	  if (tid == 0)
	      d_U[b] = sdata[0] / d;
	}
}










